This notebook runs a longitudinal-functional analysis on simulated data.
library(tidyverse)
library(MASS)
library(Rcpp)
library(splines)
library(LFBayes)
library(plotly)
### Control parameters
set.seed(999)
source("/Users/johnshamshoian/Documents/R_projects/LFBayes/Example/Simulation_funcs.R")
errorvar <- .025
SS <- 20
TT <- 20
t <- seq(from = 0, to = 1, length.out = TT)
s <- seq(from = 0, to = 1, length.out = SS)
n <- 60
tt <- list()
tt[[1]] <- 1:(TT*SS)
tt <- rep(tt, n)
p1 <- 12
p2 <- 12
q1 <- 3
q2 <- 3
Bt <- bs(t, df = p1, intercept = TRUE)
Bs <- bs(s, df = p2, intercept = TRUE)
Bt1 <- bs(t, df = p1, intercept = TRUE)
Bs1 <- bs(s, df = p2, intercept = TRUE)
### Generate key model quantities
H <- GenerateH(q1, q2)
mu1 <- GenerateMu1(s,t)
Lambda <- Loading.Matern(t, p1, q1, Bt)
Gamma <- Loading.Brown.Bridge(s, p2, q2)
Cov <- kronecker(Bs%*%Gamma, Bt%*%Lambda)%*%H%*%
t(kronecker(Bs%*%Gamma, Bt%*%Lambda)) + errorvar * diag(SS * TT)
### Generate data from model
x <- mvrnorm(n, mu = as.vector(mu1), Sigma = Cov)
sx <- sd(x)
mx <- mean(x)
x <- (x-mx)/sx
Smooth_scaled_cov <- (Cov - errorvar * diag(SS * TT)) / sx^2
mu <- (mu1 - mx)/sx
y <- lapply(1:n, function(i) x[i,])
missing <- list()
for(ii in 1:n){
missing[[ii]] <- numeric(0)
}
X <- cbind(rep(1, n))
### Visualize a few trajectories
nsub <- 6
small_data <- tibble(id = numeric(),
func_time = numeric(),
long_time = numeric(),
value = numeric())
small_data <- small_data %>%
add_row(id = rep(1:nsub, each = SS * TT),
func_time = rep(rep(t, SS), nsub),
long_time = rep(rep(s, each = TT), nsub),
value = c(t(x[1:nsub,])))
id.labs <- paste("Subject", 1:nsub)
names(id.labs) <- "1":nsub
small_data %>%
filter(long_time %in% c(s[5], s[10], s[15])) %>%
ggplot(aes(func_time, value)) +
geom_line(aes(linetype = factor(long_time))) +
facet_wrap(. ~ id, labeller = labeller(id = id.labs)) +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
labs(x = "Functional time", y = "Response", linetype = "Longitudinal time") +
scale_linetype_discrete(labels = round(c(s[5], s[10], s[15]), 2))
### MCMC control parameters
iter <- 500 # Number of iterations
burnin <- 100 # Burnin iterations
thin <- 1 # Thinning for each chain
nchain <- 1 # Number of chains
q1s <- 3 # Number of latent factors for functional dimension
q2s <- 3 # Number of latent factors for longitudinal dimension
alpha <- .05 # Type 1 error
neig <- 3 # Number of eigenfunctions for inference
### Processing
MCMC <- mcmcWeakChains(y, missing, X, Bs1, Bt1,
q1s, q2s, iter, thin, burnin, nchain)
MCMC_eigen <- eigenLFChains(Bs1, Bt1, MCMC, neig,
iter, burnin, nchain, s, t, alpha)
### Prepare data for plotting
postmean <- matrix(MCMC_eigen$postmean, nrow = TT, ncol = SS)
fig <- plot_ly(z = ~postmean, x = s, y = t) %>% add_surface()
eigenfunction_tibble <- tibble(time = numeric(), type = character(),
value = numeric(), bound = character(),
number = numeric())
eigenfunction_tibble <- eigenfunction_tibble %>%
add_row(value = c(MCMC_eigen$eigvecFuncmean,
MCMC_eigen$eigvecFunclower,
MCMC_eigen$eigvecFuncupper),
number = rep(rep(1:neig, each = TT), 3),
bound = rep(c("mean", "lower", "upper"), each = neig * TT),
type = "Functional",
time = rep(t, 3 * neig)) %>%
add_row(value = c(MCMC_eigen$eigvecLongmean,
MCMC_eigen$eigvecLonglower,
MCMC_eigen$eigvecLongupper),
number = rep(rep(1:neig, each = SS), 3),
bound = rep(c("mean", "lower", "upper"), each = neig * SS),
type = "Longitudinal",
time = rep(s, 3 * neig))
### Plot mean surface
axx <- list(
title = "Longitudinal time"
)
axy <- list(
title = "Functional time"
)
axz <- list(
title = "Posterior mean"
)
fig <- plot_ly() %>%
add_surface(z = ~postmean, x = s, y = t) %>%
layout(scene = list(xaxis=axx,yaxis=axy,zaxis=axz,
aspectratio = list(x = .9, y = .8, z = 1))) %>%
hide_colorbar()
fig
### Plot eigenfunctions
options(repr.plot.width=6, repr.plot.height=3)
number.labs <- paste("Eigenfunction", 1:neig)
names(number.labs) <- c("1":neig)
eigenfunction_tibble %>%
ggplot(aes(time, value)) +
geom_line(aes(linetype = bound)) +
facet_grid(type ~ number, labeller = labeller(number = number.labs)) +
scale_linetype_manual(values=c("dashed", "solid", "dashed"))+
theme_bw() +
theme(legend.position = "none",
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
### Percent variability explained by first few eigenfunctions in functional direction
100 * rowMeans(MCMC_eigen$eigvalFunc)[p1:(p1-neig+1)]
## [1] 63.379161 24.612260 7.079978
### Percent variability explained by first few eigenfunctions in longitudinal direction
100 * rowMeans(MCMC_eigen$eigvalLong)[p2:(p2-neig+1)]
## [1] 78.402375 14.563600 2.790165